//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#ifndef __dynamics_hpp__
#define __dynamics_hpp__

#include "math3d.hpp"
#include "utils.hpp"
#include "threading.hpp"
#include "statics.hpp"

class Statics;

/** Dynamic element for fields at mid-step calculations. */
class DynamicElement
{
public:
  virtual void addFields (const vect3d& posv, prec_t t, vect3d& bfield, vect3d& efield) = 0;
  virtual void addPotential (const vect3d& posv, prec_t t, prec_t& potential) = 0;

  virtual ~DynamicElement();
};

/** Movement of particles. */
class ElectroDynamics
{
public:
  ElectroDynamics (Statics& statics, bool disableThreads = false);
  ~ElectroDynamics ();

  /** Get bfield and efield at given position.
    * @param except_electron Particle ID to ignore in the calculations
    */
  void getFields (const vect3d& posv, vect3d& bfield, vect3d& efield, size_t except_electron = (size_t)-1);

  void getFields (const vect3d& posv, prec_t t, vect3d& bfield, vect3d& efield, size_t except_electron = (size_t)-1);

  /** Get the electric potential at the given point. */
  prec_t getPotential (const vect3d& posv, size_t except_particle = (size_t)-1);

  /** Get bfield at given position. */
  void getBField (const vect3d& posv, vect3d& bfield) { vect3d tmp; getFields (posv, bfield, tmp); }
  /** Get efield at given position. */
  void getEField (const vect3d& posv, vect3d& efield) { vect3d tmp; getFields (posv, tmp, efield); }

  /** Returns true/false whether given position is in contact with any static or dynamic elements.
    * @param except Don't count given particle ID.
    */
  bool inContact (const vect3d& pos, size_t except = (size_t)-1);

  /** Step simulation one timeunit ahead.
    * In a simulation one would continuously call this.
    */
  void step (size_t count = 1);

  /** Elapsed time in seconds of the simulation. */
  prec_t get_elapsed_time ();

  /** Inject a particle into the simulation.
    * @param mass Mass of particle (kg)
    * @param charge Charge of particle (coulomb)
    * @param posv Position of injection
    * @param heading Particle heading
    * @param electronVolts Particle speed/energy
    * @return Particle ID.
    */
  size_t inject_particle (prec_t mass, prec_t charge, const vect3d& pos, const vect3d& heading, prec_t electronVolts);

  /** Inject an electron into the simulation.
    * @param posv Position of injection
    * @param heading Particle heading
    * @param electronVolts Particle speed/energy
    * @return Particle ID.
    */
  size_t inject_electron (const vect3d& posv, const vect3d& heading, prec_t electronVolts);

  /** Update particle's position and speed. */
  void update (size_t index, const vect3d& posv, const vect3d& heading, prec_t electronVolts);

  /** Reset particles' current field condition. */
  void reset ();

  /** Erase a particle with given index. */
  void erase (size_t index);

  // electric field force: qE
  // lorentz force q * v x B
  // total force: q * (E + v x B)
  // electric field module of a charge: (1/(4*pi*e0)) * (Q/|r|*|r|) * r^
  // electric field: rvect * Ke / (rvect*rvect)^(3/2)
  // bfield of a single electron: B = Km * q * (v x r) / (r * r)^(3/2)
  // F = ma; a = F/m
  // delta_v = a * delta_t = F * delta_t / m

  /** Debug all particles, their positions, speed, etc. */
  void print ();

  /** Get position of the given particle. */
  const vect3d& get_pos (size_t i) const { return particles[i].pos; }
  /** Get speed (m/s) of the given particle. */
  const vect3d& get_v (size_t i) const { return particles[i].v; }
  /** Get charge of the given particle. */
  const prec_t get_q (size_t i) const { return particles[i].charge; }

  /** Set simulation accuracy and timestep.
    * @param accuracy Substep adaptation level, bigger means adapt more means slower simulation
    * @param stepf Step size modifier, bigger means simulation of smaller timesteps - need to call more step()-s to simulate the same time
    */
  void setDTM (prec_t accuracy, prec_t stepf = 1.0);
  /** @deprecated */
  void setPrintStream (FILE *fp) { printfp = fp; }

  /** Convert electronVolts into speed (meters per second) for given electronVolts and mass. */
  static prec_t eV2ms (prec_t ev, prec_t mass);
  /** Convert speed (meters per second) into electronVolts for given speed and mass. */
  static prec_t ms2eV (prec_t ms, prec_t mass);

  /** Count of particles in the system. */
  size_t size () { return particles.size(); }

  Statics& get_statics () { return statics; }

  /** How many steps have passed since beginning. */
  size_t get_ticks () { return ticks; }

  /** How many ticks of simulation would take for given speed (meters per second) to cover the distance. */
  size_t ticks_for_distance (prec_t speed, prec_t distance);

  /** Disable electron-electron fields calculations. */
  void setDisableFields (bool df) { disable_fields = df; }

  /** Save dynamics state into given file. */
  bool save (const char *path);

  /** Load dynamics state from the given file. */
  bool load (const char *path);

  void addDynamicElement (DynamicElement *elem);

private:
  Statics& statics;
  size_t ticks;
  prec_t delta_t;
  prec_t adaptf, madaptf;
  FILE *printfp;

  bool disable_fields;

  struct particle
  {
    prec_t mass;
    prec_t charge;
    prec_t charge_by_mass;
    vect3d pos;
    vect3d v;
    vect3d next_pos;
    vect3d next_v;
    vect3d last_bf;
    vect3d last_ef;
    size_t adaptlevel;
  };
  typedef std::vector<particle> particles_t;
  particles_t particles;

  struct pp
  {
    size_t count;
    prec_t dt;
    prec_t hdt;
    prec_t dt_by_6;
    prec_t dtdt;
    prec_t dtdt_by_2;
    prec_t dtdt_by_6;
    prec_t dtdt_by_8;

    prec_t qdt;
    prec_t qhdt;
  };
  pp *pmap;

#ifdef HAVE_THREADS
  TaskManager *tmgr;
  struct steptask : public Task
  {
    ElectroDynamics *dyn;
    size_t n, c;
    void execute () { dyn->run_steptask (n, c); }
    virtual ~steptask () { }
  };
  friend struct steptask;

  void run_steptask (size_t n, size_t c);

  std::vector<steptask> tasks;
#endif

  void internal_step (size_t i);
  void dynstep (size_t i);

  void update_pmap ();

  typedef std::vector<DynamicElement*> delems_t;
  delems_t delems;
};

/** Dynamic element consisting of three static elements extrapolating fields
  * using 3-point Lagrange polynomial interpolation. */
class Dynamic3 : public DynamicElement
{
public:
  Dynamic3 ();
  virtual ~Dynamic3();

  void next ();
  Statics& get(int n) { int i = index + n; if ( i > 2 ) i -= 3; return statics[i]; }

  int head() const { int i = index + 2; if ( i > 2 ) i -= 3; return i; }

  void addFields (const vect3d& posv, prec_t t, vect3d& bfield, vect3d& efield);
  void addPotential (const vect3d& posv, prec_t t, prec_t& potential);

public:
  int index;
  Statics statics[3];

  Statics *s0, *s1, *s2;
};

#endif // !__dynamics_hpp__

